data {
  m.1 <- m - 1
}
model {
  for (i in 1:n) {
    y[i, 1:m] ~ dmulti(theta[i, 1:m], k[i])
    for (j in 1:m.1) {
      theta[i, j] <- exp.theta.star[i, j] / sum.exp.theta.star[i]
    }
    theta[i, m] <- 1 / sum.exp.theta.star[i]
    sum.exp.theta.star[i] <- sum(exp.theta.star[i, 1:m.1]) + 1
    for (j in 1:m.1) {
      exp.theta.star[i, j] <- exp(theta.star[i, j])
    }
  }
  for (i in 2:n.1) {
    for (j in 1:m.1) {
      theta.star[i, j] ~ dnorm(mu[i, j], Omega.star[j])
      mu[i, j] <- theta.star[i - 1, j] + inprod(beta[j, ], x[i, ])
    }
  }
  theta.star[1, 1:m.1] ~ dmnorm(mu0, Omega.star0)
  for (i in (n.1 + 2):n) {
    for (j in 1:m.1) {
      theta.star[i, j] ~ dnorm(mu[i, j], Omega.star[j])
      mu[i, j] <- theta.star[i - 1, j] + inprod(beta[j, ], x[i, ])
    }
  }
  theta.star[n.1 + 1, 1:m.1] ~ dmnorm(mu0, Omega.star0)
  for (i in 1:m.1) {
    for (j in 1:l) {
      beta[i, j] ~ dnorm(0, 0.01)
    }
    Omega.star[i] ~ dscaled.gamma(100, 2)
    Omega[i] <- pow(Omega.star[i], -2)
  }
}